26/06, 2020

Estudio de abundancia

library(rworldxtra)
library(raster)
library(sf)
library(tidyverse)

Datos son parte de PREDICTS

data("countriesHigh")
Datos <- read_csv("https://raw.githubusercontent.com/derek-corcoran-barrios/derek-corcoran-barrios.github.io/master/Presentaciones_Espacial/Bombus.csv")
## Miremos los datos
Datos <- Datos %>% st_as_sf(coords = c(5, 6), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
Mapa <- countriesHigh %>% st_as_sf() %>% st_crop(Datos)
write_sf(Datos, "Datos.shp")

Vamos graficando

ggplot() + geom_sf(data = Mapa) + theme_bw()

Agreguemos los puntos

ggplot() + geom_sf(data = Mapa) + geom_sf(data = Datos) + theme_bw()

Agreguemos colores por especie

ggplot() + geom_sf(data = Mapa) + geom_sf(data = Datos, aes(color = species)) + 
    theme_bw()

Agreguemos el tamaño por abundancia

ggplot() + geom_sf(data = Mapa) + geom_sf(data = Datos, aes(color = species, 
    size = Measurement)) + theme_bw()

o un panel por especie

ggplot() + geom_sf(data = Mapa) + geom_sf(data = Datos) + theme_bw() + 
    facet_wrap(~species)

Resumen

species n mean max min
Bombus impatiens 400 7.8675000 106 0
Bombus bifarius 382 7.2251309 131 0
Bombus bimaculatus 400 2.5850000 52 0
Bombus vosnesenskii 382 2.3612565 81 0
Bombus pensylvanicus 382 1.3926702 67 0
Bombus griseocollis 18 0.5555556 6 0
Bombus auricomus 18 0.4444444 8 0
Bombus occidentalis 382 0.3376963 23 0
Bombus fervidus 18 0.2777778 2 0
Bombus terricola 382 0.0811518 13 0
Bombus affinis 382 0.0575916 12 0

Ahora a modelar

Modelemos Bombus impatiens

B_impatiens <- Datos %>% dplyr::filter(species == "Bombus impatiens")
ggplot() + geom_sf(data = Mapa) + geom_sf(data = B_impatiens, 
    aes(size = Measurement)) + theme_bw()

Obtengamos datos climaticos

Bioclim <- getData("worldclim", var = "bio", res = 2.5) %>% crop(B_impatiens)
plot(Bioclim)

Que es esto?

Ahora juntemos las bases de datos

Extracción de datos en los puntos

Clima <- extract(Bioclim, B_impatiens) %>% as.data.frame()
bio1 bio7 bio12 bio15
62 396 507 65
161 378 896 39
118 405 1004 33
107 416 934 39
118 392 958 22
118 392 961 22
109 385 970 24
75 436 849 49
97 419 904 41
76 435 867 44
124 356 1243 13
72 400 581 56
68 406 476 68
47 475 413 67
48 477 440 63
56 486 424 65
125 435 602 54
160 386 755 44
196 320 833 32
127 415 816 42
200 324 850 32
168 359 1035 30
124 401 923 47
208 258 1201 30
68 453 739 51
137 353 1218 22
128 379 1110 25
182 318 1187 15
129 379 1098 25
119 397 983 28
30 460 712 48
73 444 821 46
66 447 826 45
110 406 943 26
109 403 936 26
109 403 936 26
109 403 936 26
69 435 815 42
69 435 815 42
77 420 814 39
131 380 1126 17
87 409 910 35
196 276 1608 18
111 393 986 23
111 393 986 23
110 393 984 24
110 393 984 24
135 366 1214 15
170 325 1460 19
104 373 958 20
133 346 1173 11
75 366 1220 17
71 376 1191 15
71 376 1201 16
75 362 1246 18
64 386 1176 15
83 403 889 32
77 420 814 39
173 358 729 36
89 436 612 55
93 440 636 58
199 323 858 32
160 373 892 40
104 428 736 53
121 414 864 48
125 407 925 45
88 428 767 52
126 396 966 44
208 258 1201 30
187 315 1291 16
143 373 1159 25
67 456 746 51
177 333 1376 16
96 421 915 37
107 406 936 38
123 395 995 32
71 455 782 48
67 442 838 46
166 343 1268 20
75 441 808 45
68 448 821 45
125 397 990 19
125 397 990 19
123 398 965 21
160 352 1309 21
85 417 878 39
65 450 830 43
128 387 983 19
118 389 954 22
69 435 815 42
98 405 935 30
101 404 930 30
77 420 814 39
131 380 1126 17
87 409 910 35
118 392 958 22
134 369 1193 15
140 364 1257 14
110 393 984 24
168 327 1419 21
114 392 994 22
111 389 986 23
89 404 918 32
198 264 1530 21
142 365 1366 15
141 360 1260 14
105 382 1045 22
106 382 1051 21
181 310 1319 22
98 391 955 25
142 357 1311 15
122 368 1170 17
179 314 1441 20
117 369 1102 17
116 369 1107 17
145 340 1432 16
125 359 1161 16
150 338 1407 15
131 338 1513 14
98 386 883 21
148 325 1457 16
133 334 1427 14
161 326 1278 15
99 375 858 22
118 341 1278 13
130 339 1137 15
126 340 1151 13
164 332 1190 15
84 302 1681 9
92 363 957 20
117 329 1378 9
114 327 1396 9
158 337 1224 13
99 364 1025 18
108 334 1164 12
168 319 1135 21
133 342 1219 10
162 331 1178 18
160 336 1173 16
130 345 1037 11
147 339 1142 12
86 359 864 13
131 353 1107 9
78 370 1043 18
95 361 987 15
85 367 845 22
83 369 880 22
53 397 1075 14
98 377 1086 10
79 389 1146 9
109 329 1149 8
92 376 1221 7
58 389 1307 8
108 329 1145 9
87 394 1146 7
60 408 1142 10
92 374 1169 7
74 384 1153 7
94 371 1150 9
85 376 865 30
97 423 889 43
82 463 451 63
41 491 467 63
69 463 584 55
79 461 616 54
66 454 665 50
68 457 777 49
112 404 976 36
125 397 990 19
89 404 918 32
89 404 918 32
105 382 1045 22
122 367 1170 17
122 367 1171 17
98 391 958 24
117 369 1102 17
125 397 990 19
104 428 734 52
102 427 800 51
97 278 1678 69
17 410 583 20
51 305 914 63
51 305 914 63
51 305 914 63
1 419 545 16
1 419 545 16
72 367 789 71
76 333 864 65
76 333 864 65
NA NA NA NA
57 316 1035 62
-1 359 566 27
47 308 1202 62
97 359 509 26
-14 343 631 25
44 354 613 38
82 370 742 72
11 399 664 18
38 382 340 15
51 371 451 28
51 371 451 28
-9 359 638 11
45 257 1990 61
43 388 682 37
-16 378 662 16
57 320 812 61
17 400 448 17
76 333 864 65
64 353 595 58
54 346 608 29
52 345 621 31
28 412 597 21
76 361 353 36
35 372 457 30
99 189 615 49
41 382 379 27
86 255 2130 62
NA NA NA NA
77 218 1466 62
62 296 644 69
9 365 518 30
86 337 697 65
95 376 407 53
96 237 1563 47
28 291 1371 58
29 321 895 58
56 312 1178 66
51 416 335 27
16 419 645 21
16 342 932 34
39 335 761 43
27 385 513 17
27 385 513 17
71 292 625 71
74 296 610 72
80 299 587 72
57 381 454 46
60 378 394 26
39 371 515 31
24 381 481 36
72 313 1331 64
20 364 498 24
18 418 355 47
38 376 399 27
98 191 669 49
93 240 1740 62
50 317 1195 67
51 359 411 43
31 324 965 56
19 363 588 18
41 382 380 28
37 366 683 20
33 413 317 31
82 422 322 62
91 260 1852 57
15 308 1028 52
63 318 865 63
16 342 932 34
25 350 900 34
16 342 932 34
16 342 932 34
40 340 876 39
4 414 625 17
32 420 586 27
6 312 822 44
83 345 437 39
5 413 669 21
26 348 520 23
35 353 477 25
51 366 397 26
51 366 397 26
35 353 477 25
52 341 733 43
31 374 512 24
31 374 512 24
31 374 512 24
9 360 591 21
16 365 563 22
97 192 603 45
48 362 726 22
52 377 567 50
69 256 2060 54
65 281 1181 70
65 281 1181 70
28 396 445 25
18 388 490 26
84 319 529 70
55 321 807 63
55 321 807 63
112 307 1295 79
77 361 396 41
41 333 778 45
40 335 773 44
98 368 646 63
61 325 1135 66
100 241 1643 65
94 189 695 55
-12 347 622 25
39 372 445 23
103 249 1322 63
11 394 441 24
38 362 444 36
52 369 591 26
96 374 404 54
60 389 460 59
19 381 720 33
72 313 1331 64
33 359 427 21
-8 351 607 26
35 385 526 29
24 389 560 24
45 257 1990 61
2 362 557 27
100 198 2482 62
22 362 570 20
43 399 439 51
70 406 439 46
58 398 361 18
28 412 597 21
49 406 362 27
25 350 900 34
46 398 456 55
58 392 308 32
45 311 1113 65
18 392 381 34
41 382 380 28
11 368 511 31
60 395 421 29
48 386 659 26
84 329 535 37
65 243 2155 52
15 314 964 52
44 375 485 30
-6 409 628 17
80 315 550 71
45 400 450 55
9 360 591 21
2 339 568 13
64 356 484 69
40 392 451 15
102 255 1481 58
17 355 608 38
24 336 762 51
25 405 339 26
55 251 2146 61
91 236 2185 64
92 236 2152 64
83 328 539 38
83 328 539 38
12 345 513 17
36 415 420 11
10 363 526 32
6 354 510 8
35 266 1918 63
27 434 482 21
61 473 308 46
20 364 498 24
73 444 348 63
85 227 2363 67
32 406 451 11
12 375 581 17
12 375 581 17
114 350 986 78
42 391 400 26
93 195 764 50
66 388 346 56
37 355 667 37
42 362 634 37
42 362 634 37
42 359 637 37
43 360 572 27
43 360 572 27
27 362 602 29
72 292 617 71
41 211 2400 59
109 342 397 33
91 209 851 62
97 192 603 45
98 368 646 63
94 189 695 55
35 258 2052 62
-17 337 736 15
99 379 918 27
99 378 913 26
101 381 931 26
93 379 902 27
96 376 902 27
99 379 918 27
99 378 913 26
99 379 917 27
99 381 917 27
101 378 918 26
100 378 913 26
98 377 907 27
89 382 891 29
85 380 872 30
97 402 929 29
96 384 916 28
93 379 902 27
92 378 899 28

Juntamos con abundancia

B_impatiens <- B_impatiens %>% bind_cols(Clima)

Veamos relaciones

A modelar

GLM

  • Para mas información ver este video
Fit1 <- glm(Measurement ~ bio1 + I(bio1^2) + bio7 + I(bio7^2) + 
    bio12 + I(bio12^2) + bio15 + I(bio15^2), family = poisson, 
    data = B_impatiens)
term estimate std.error statistic p.value
(Intercept) 3.0529399 1.1192789 2.7275954 0.0063798
bio1 0.0559696 0.0035197 15.9017989 0.0000000
I(bio1^2) -0.0002341 0.0000154 -15.1665679 0.0000000
bio7 -0.0731175 0.0065246 -11.2063574 0.0000000
I(bio7^2) 0.0001204 0.0000094 12.8028406 0.0000000
bio12 0.0126855 0.0008787 14.4364480 0.0000000
I(bio12^2) -0.0000048 0.0000004 -12.5696430 0.0000000
bio15 -0.0409474 0.0070566 -5.8027476 0.0000000
I(bio15^2) -0.0000698 0.0001102 -0.6334361 0.5264489

Predecir en un raster

Prediccion <- predict(Bioclim, Fit1, type = "response")
plot(Prediccion, colNA = "black")

Dudas?

Mapa con SF

Prediccion_DF <- Prediccion %>% as("SpatialPixelsDataFrame") %>% 
    as.data.frame() %>% rename(Abundancia = layer)
ggplot() + geom_tile(data = Prediccion_DF, aes(x = x, y = y, 
    fill = Abundancia)) + geom_sf(data = Mapa, alpha = 0) + theme_bw() + 
    scale_fill_viridis_c() + xlab("") + ylab("")

Mapa con los puntos

ggplot() + geom_tile(data = Prediccion_DF, aes(x = x, y = y, 
    fill = Abundancia)) + geom_sf(data = Mapa, alpha = 0) + geom_sf(data = B_impatiens, 
    aes(size = Measurement)) + theme_bw() + scale_fill_viridis_c(option = "plasma") + 
    theme(legend.position = "bottom") + xlab("") + ylab("")

Mapa con los puntos

Predicción futura

Capas

Futuro <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, 
    model = "HD", year = 70) %>% crop(B_impatiens)
Futuro <- Futuro[[c(1, 7, 12, 15)]]
names(Futuro)
## [1] "hd85bi701"  "hd85bi707"  "hd85bi7012" "hd85bi7015"
names(Futuro) <- c("bio1", "bio7", "bio12", "bio15")
Futuro_pred <- predict(Futuro, Fit1, type = "response")
plot(Futuro_pred, colNA = "black")

Capas

Comparacion